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We discuss numerical complexity of the Liischer algorithm applied to the Hubbard Model. In particular we 
present comparison to a certain algorithm, based on direct computation of the fermionic determinant. 
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1. Introduction 

Recently a new type of algorithms for nu- 
merical simulations of models with dynamical 
fermions has been proposed by M. Liischer [Q. 
The idea is based on the approximation of a 
fermionic determinant by the path integral over 
some number of bosonic fields, which allows ones 
to elegantly avoid the main difficulty coming from 
the nonlocality of the fermionic determinant. The 
complexity of a sweep of the Liischer algorithm 
is proportional to the volume of the system and 
in this respect it wins in comparison with the 
standard hybrid Monte Carlo Of course, 

what is interesting in practice is a comparison of 
the computer cost needed to produce indepen- 
dent configurations which is the product of the 
sweep cost and the autocorrelation time. There- 
fore the crucial question is how the autocorre- 
lation time depends on the volume of the sys- 
tem. Contrary to the sweep time the autocorre- 
lation time can strongly depend on the external 
parameters. Therefore it is particularly interest- 
ing to test the dependence of the autocorrelation 
time on the various interesting physical parame- 
ters. Recent results obtained in QCD are rather 
promising showing that the computational cost 
needed to obtain independent configurations is 
comparable to the standard HMC algorithms || . 

In this letter we present new results of numer- 
ical simulations of the Hubbard Model using the 
Liischer method and compare them to standard 
methods in context of the algorithmical complex- 
ity. This is a continuation of our contribution 
to the proceedings Lattice 95 where we presented 
preliminary results H. In particular, we noted 



there possible improvements which reduce auto- 
correlation times. Here, we study the errors intro- 
duced by the bosonic approximation and compare 
dynamical properties of the multiboson algorithm 
with the simple exact algorithm. 

2. Liischer approximation 

We consider the Hubbard Model defined by the 
Hamiltonian 



H = 



K a ™ a i° 



(i) 



where fermionic operators acting on elec- 

trons on site i of the lattice. The sum runs 
over nearest neighbors (each symmetric pair (i, j) 
is counted twice). Physical parameters entering 
the model arc: K - hopping parameter, U - the 
strength of the effective Coulomb interaction, and 
the inverse temperature (3. The interaction term 
of the Hamiltonian describes the system with the 
half-filled band. One find the following represen- 
tation of the partition function for the model M : 



Z ~ / [dAd(j>] exp 
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where Q)Q = M^M/X ma x and X ma x denotes the 
largest eigenvalue of the matrix M^Ai. The ma- 
trix M. is the fermion matrix entering the Eu- 
clidean formulation 
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A(x) is the Hubbard-Stratanovich continuous 
field, and <j>k(x) are the auxiliary bosonic fields. 
The particular form of the polynomial P n {x) 
used in the approximation of the function 1/x, is 
built from the Chebyshev polynomials and gives 
a rapid uniform convergence to the function 1/x 
in the interval (e, 1). Its roots are introduced al- 
ready in (||) and their real and imaginary parts 
are denoted by + ifiki k = 1, n respectively. 
The relative error of the approximation defined as 
|-R(a:)| = \xP n {x) — 1| is exponentially bounded 

|B(I)I<2 (^)" +1 ' (3) 

One sees that to decrease e one has to increase 
n. The most economic choice of e corresponds 
to the largest possible value which still gives a 
good approximation of the determinant. To make 
this statement precise we introduce following the 
quantity H as a measure of error (or goodness) 



8 = \y 1/v 

where 
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det (Q t g)P n (Q t Q). 



(4) 



(5) 



The power 1/V, where V is naturally adjusted to 
the dimension of the matrix M.. Figure 1 shows 
5 as a function of e for a different numbers of 
fields measured on a one typical configuration. 
The optimal value of e lies close to the smallest 
eigenvalue of the matrix Q^Q (X m in — 0.032). 
Such a behavior is consistent with the previous 
findings in the QCD case. 

3. Exact algorithm 

To compare the computational cost of the 
multiboson algorithm with the cost of the more 
standard approaches, we implemented a simple 
exact algorithm based on the updating scheme al- 
lowing for computing changes in the matrix M~ x 
while making trial changes in the matrix M.. In 
order to update the field A located on the site i 
one must calculate the ratio of the fcrmion deter- 
minants. Since the only A dependent elements of 




Figure 1. Errors of the Luscher approximation 
measured on a single typical configuration. The 
quantity 5 defined in text is shown as a function of 
e. The solid, dashed and dotted line is for number 
of fields 50, 70 and 100 respectively 



the matrix M lie on the diagonal we consider the 
following change in the matrix M. , 



M' = (I + A)M, 



(G) 



where A is matrix with one nonzero element 



dijd. Then 



dctTW' 
dctTW 



= dot (/ + M~ 1 A) = 1 + M^d 



(7) 



is determined completely by elements of the ma- 
trix Mr 1 . Once the trial change has been ac- 
cepted the updated matrix Mr 1 can be evaluated 
from the formula 



M 



'-i 



1 + dMZ 



(8) 



This process is very economical from the com- 
putational point of view since one update of M _1 
requires 0{V 2 ) operations comparing with 0(V 3 ) 
operations needed to evaluate the determinant 
with the brute force method. Also the today's 
leading algorithm for simulations of the Hubbard 
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Figure 2. Eigenvalues of the matrix M. on the 
complex plane 



Model is essentially based on it However, con- 
trary to the standard one our particular formula- 
tion of path integrals does not require multiplica- 
tion badly conditioned matrices. Thus we believe 
that this is free of the numerical instabilities. In- 
deed, we performed thousands sweeps at (3 up 
to 8 without accumulating numerical errors. An 
additional cost comes from working with larger 
matrices. 

4. Final remarks 

As we have previously reported the strong au- 
tocorrelation between generated configurations is 
a serious obstacle which limits the available range 
of parameters to U < 2 and (3 = 1 for the multi- 
boson algorithm. In contrast the program based 
on exact evaluation of the fcrmionic determinant 
has no such restrictions and works equally good 
at U = 1 and U = 4. 

The comparison of CPU time needed for both 
algorithms varies accordingly to the lattice sizes. 
On small lattices exact algorithm is substantially 
faster. For example for a lattice 6 2 8 in our exper- 
iments on the workstation HP735/125 it takes 20 
minutes which should be compared with 10 hours 



needed for our implementation of the Liischer al- 
gorithm to produce results of the same quality. Of 
course for larger lattices the computational cost 
grows much more faster for the exact algorithm. 
In the region of weak coupling we were able to 
could simulate with the Liischer algorithm lat- 
tices with the size up to 16 2 8. 

The polynomial approximation of the function 
1/x can be extended to the complex plane ||. 
Thus one try to approximate directly the inverse 
of the matrix M. This would reduce the condi- 
tion numbers for matrices entering the problem 
and would result in the simpler final action. This 
modification of the original Liischer idea gave 
very promising results However, in contrast 
to the QCD case the matrix M. for the Hubbard 
Model has eigenvalues with the positive and neg- 
ative real parts as can for example be seen in Fig- 
ure 2. Because one cannot extend the domain of 
the convergence beyond the singular point (0,0), 
this is, unfortunately, impossible to adopt this 
modification to the Hubbard Model. 
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